library(ggplot2)
library(stargazer)
library(MASS)



## 1: Analyses related to the influx of outsiders ======

# Get the base data
crossmove_data <- read.csv("crossmove_data.csv", header = T)

# Note: For columns Position_From and Position_To, 1 = Periphery and 0 = Center

# == == == == === == === == == == == == === === == ==  == == == == == === == == ==
# == == == ==                                                              = == ==
# == == == == Table 1: Movements Involving the Center (C)                  == == = 
# == == == ==          and Periphery (P)                                   = == ==
# == == == ==                                                              = == ==
# == == == == === == === == == == == == === === == ==  == == == == == === == == ==

# Table enumerating movements involving center and periphery
# Periphery = 1 and Center = 0
temp <- table(crossmove_data$Position_From, crossmove_data$Position_To)
moves <- c(temp[1], temp[3], temp[2], temp[4])

# Table enumerating movements involving military, parliament, etc.
temp_type <- table(crossmove_data$Position_From_Name, crossmove_data$Position_To_Name)

type_sum <- c( sum(temp_type[2, ], temp_type[ , 2]),
               sum(temp_type[5, ], temp_type[1:4, 5]),
               sum(temp_type[4, ], temp_type[ , 4]),
               sum(temp_type[3, ], temp_type[ , 3]))


# Table 1
table_1 <- cbind(data.frame(Type_of_Movement = c("C<->C", "C->P", "P->C", "P<->P"),
                            Amount = moves,
                            Proportion = moves/nrow(crossmove_data)),
                 data.frame(Periphery_Type = c("Embassy", "Province", "Parliament", "Military"),
                            Amount = type_sum,
                            Proportion = type_sum/nrow(crossmove_data)))


table_1



rm(table_1, moves, temp_type, type_sum, temp)


# == == == == === == === == == == == == === === == ==  == == == == == === == == ==
# == == == ==                                                         === == == ==
# == == == == Figure 1: Circulation of officials involving the center == == ==  == 
# == == == ==                                                         === == == ==
# == == == == === == === == == == == == === === == ==  == == == == == === == == ==


# Number of appointments broken down by Periphery-Center-Year
periph_min <- as.data.frame(table(crossmove_data[which(crossmove_data$Position_From == 1 &
                                                        crossmove_data$Position_To == 0) , 8]))

min_min <- as.data.frame(table(crossmove_data[ which(crossmove_data$Position_From == 0 &
                                                      crossmove_data$Position_To == 0) , 8 ]))

min_periph <- as.data.frame(table(crossmove_data[which(crossmove_data$Position_From == 0 &
                                                         crossmove_data$Position_To == 1) , 8]))

periph_periph <- as.data.frame(table(crossmove_data[which(crossmove_data$Position_From == 1 &
                                                           crossmove_data$Position_To == 1) , 8]))

# Construct data frame (i.e., move_descript) where data will be collated
base <- rep(0, length(42:74))
move_descriptive <- as.data.frame(cbind(42:74, base, base, base))
rm(base)

# Populate move_descript with values for relevant years and movement 
move_descriptive[which(move_descriptive$V1 %in% min_min$Var1), 2] <- min_min$Freq
move_descriptive[which(move_descriptive$V1 %in% min_periph$Var1), 3] <- min_periph$Freq
move_descriptive[which(move_descriptive$V1 %in% periph_min$Var1), 4] <- periph_min$Freq

colnames(move_descriptive) <- c("year", "min_min", "min_periph", "periph_min")

rm(periph_min, periph_periph, min_periph, min_min)
  
# == == == === == == == Plot the Results == == === == === == == === == === == #


ggplot(move_descriptive, aes(x = year, y = rowSums(move_descriptive[ , 2:3]))) +
  geom_point(color = "lightgray") +
  geom_line(color = "lightgray", linetype = "dashed") +
  geom_point(aes(x = year, y = periph_min)) +
  geom_line(aes(x = year, y = periph_min)) +
  scale_x_continuous(breaks = pretty(move_descriptive$year, n = 8)) +
  geom_vline(xintercept = 60, linetype = "dashed", size = 0.2) +
  annotate("text", 
           x = 57.5, y = 47, 
           size = 1.8,
           label = "Failed \n Coup") +
  geom_segment(aes(x = 58.5, y = 47, xend = 59.5, yend = 47), 
               size=0.2, 
               arrow = arrow(length = unit(0.02, "npc"))) +
  annotate("text", 
           x = 64, y = 27,
           size = 1.8, 
           label = "Periphery to \n Center") +
  geom_curve(aes(x = 63, y = 24.5, xend = 62, yend = 19), 
             size=0.2, 
             curvature = -0.2,
             arrow = arrow(length = unit(0.02, "npc"))) +
  annotate("text", 
           x = 67, y = 45,
           size = 1.8, 
           label = "Center to Center \n + \n Center to Periphery") +
  geom_curve(aes(x = 68, y = 40, xend = 70, yend = 35), 
             size=0.2, 
             curvature = -0.2,
             arrow = arrow(length = unit(0.02, "npc"))) +
  labs(y = "Number of Appointments", x = "Year") +
  theme_classic()

rm(move_descriptive)
# == == == == === == === == == == == == === === == ==  == == == == == === == == ==
# == == == ==                                                        = === == == =
# == == == == Table 2: Determinants for Periphery to Center Movement == == == == =
# == == == ==                                                        = === == == =
# == == == == === == === == == == == == === === == ==  == == == == == === == == ==


# Prune dataset to only movements emanating from the periphery
periph <- crossmove_data[which(crossmove_data$Position_From == 1), ]

# Run regression in accordance with equation 1 in the paper.
# Each model contains more control variables.
periph_logit.1 <- glm(move_center ~ coup, 
                      data = periph, family = "binomial")

periph_logit.2 <- glm(move_center ~ coup + 
                        Return + 
                        Seq + 
                        Time_Past, 
                      data = periph, family = "binomial")


periph_logit.3 <- glm(move_center ~ coup + 
                        Return + 
                        Seq + 
                        Time_Past + 
                        embassy + province + parliament + military,
                      data = periph, family = "binomial")

periph_logit.4 <- glm(move_center ~ coup + 
                        Return + 
                        Seq + 
                        Time_Past + 
                        embassy + province + parliament + military +
                        coup*Return, 
                      data = periph, family = "binomial")


# Output table
stargazer(periph_logit.1, periph_logit.2, periph_logit.3, periph_logit.4, 
          covariate.labels = c("Coup", "Returnee", "Sequence", "Time Elapsed",
                               "Embassy", "Province", "Parliament", "Military",
                               "Coup:Returnee", "Constant"),
          style = "AJPS", 
          type = "text")

# clear 
rm(periph_logit.1, periph_logit.2, periph_logit.4)


# == == == == === == === == == == == == === === == ==  == == == == == === == == ==
# == == == == === == === == ==                               == == == === == == ==
# == == == == === == === == ==  Figure 2: First difference  == == == == == == == == 
# == == == == === == === == ==                               == == == === == == ==
# == == == == === == === == == == == == === === == ==  == == == == == === == == ==

# Using model periph_logit.3, simulate draws from the multivariate distribution
# to obtain coefficients with which to calculate the first difference and 
# associated confidence interval

set.seed(231)

# 1000 draws from the the multivariate distribution (i.e., variance-covariance matrix)
sim_logit_coef <- mvrnorm(1000, coef(periph_logit.3), vcov(periph_logit.3))

# Set up new data frame (i.e., dv_set) with values set at their mean or median
dv_set <- cbind(Intercept = 1, 
                with(periph, 
                     data.frame(coup = 0, Return = 0, 
                                Seq = mean(Seq), Time_Past = median(Time_Past),
                                embassy = 0, province = 0,
                                parliament = 0, military = 0)))

dv_set <- rbind(dv_set, dv_set, dv_set, dv_set, dv_set, dv_set)

# For each row, turn the relevant variable on
dv_set[2,2] <- 1 # Coup = 1
dv_set[3,6] <- 1 # Embassy = 1
dv_set[4,7] <- 1 # Province = 1
dv_set[5,8] <- 1 # Parliament = 1
dv_set[6,9] <- 1 # Military = 1

# Combine with the simulated coeficcients 
reference <- sim_logit_coef %*% t(dv_set[1, ])
coup_1 <- sim_logit_coef %*% t(dv_set[2, ])
embassy_1 <- sim_logit_coef %*% t(dv_set[3, ])
province_1 <- sim_logit_coef %*% t(dv_set[4, ])
parliament_1 <- sim_logit_coef %*% t(dv_set[5, ])
military_1 <- sim_logit_coef %*% t(dv_set[6, ])

# Convert Logg odds into probabilities
coup_diff <- (exp(coup_1) / (1+ exp(coup_1))) - (exp(reference) / (1+ exp(reference)))
embassy_diff <- (exp(embassy_1) / (1+ exp(embassy_1))) - (exp(reference) / (1+ exp(reference)))
province_diff <- (exp(province_1) / (1+ exp(province_1))) - (exp(reference) / (1+ exp(reference)))
parliament_diff <- (exp(parliament_1) / (1+ exp(parliament_1))) - (exp(reference) / (1+ exp(reference)))
military_diff <- (exp(military_1) / (1+ exp(military_1))) - (exp(reference) / (1+ exp(reference)))

# Prepare predicted probabilities for visualization
diff_data <- as.data.frame(cbind(rbind(mean(coup_diff), mean(embassy_diff), 
                                       mean(province_diff), mean(parliament_diff), 
                                       mean(military_diff)), 
                                 rbind(1.96 * sd(coup_diff), 1.96 * sd(embassy_diff), 
                                       1.96 * sd(province_diff),1.96 * sd(parliament_diff), 
                                       1.96 * sd(military_diff))))

diff_data <- cbind(variable = c("Coup", "Embassy", "Province", "Parliament", "Military" ),
                   diff_data)

colnames(diff_data) <- c("Variable", "Mean", "CI")

rm(reference, coup_1, embassy_1, province_1, parliament_1, military_1,
   coup_diff, embassy_diff, province_diff, parliament_diff, military_diff)

# Plot the graph
ggplot(diff_data[-5, ], aes(x = Variable, y = Mean)) +
  geom_pointrange(aes(ymin = Mean - CI, ymax = Mean + CI)) +
  geom_point(size = 2, fill = "black",  shape = 21) +
  geom_hline(yintercept = 0, linetype = 2, color = "grey60") +
  ylab("Difference in Predicted Probability") + xlab("") +
  theme_classic() 

# clear
rm(diff_data, dv_set, periph, periph_logit.3, sim_logit_coef)

# == == == == === == === == == == == == === === == ==  == == == == == === == == ==
# == == == == === Figure 3: Logit regression with the probability   == == == == == 
# == == == == ===           of crosscutting movements (i.e. C<->P)  == == == == ==
# == == == == === == === == == == == == === === == ==  == == == == == === == == ==


ratio_data <- crossmove_data

# Run the logit model
logit_ratio <- glm(cross_move ~ c_official + coup + 
                     c_official*coup + Seq + Time_Past,
                   data = ratio_data, 
                   family = "binomial")

# Create data frame to visualize coefficients from logit_ratio model
temp <- c("Intercept","C. Official","Coup", "Sequence",
          "Time Elapsed", "C. Official*Coup")

coef_plot_data <- as.data.frame(cbind(temp, 
                                      summary(logit_ratio)$coefficients[ , -c(3,4)], 
                                      confint(logit_ratio)))

colnames(coef_plot_data) <- c("Names","Estimate", "Std. Error", "Low", "High")
row.names(coef_plot_data) <- c()
coef_plot_data$Estimate <- as.numeric(as.character(coef_plot_data$Estimate))
coef_plot_data$Low <- as.numeric(as.character(coef_plot_data$Low))
coef_plot_data$High <- as.numeric(as.character(coef_plot_data$High))


# Plot graph 
ggplot(coef_plot_data[-1, ], aes(x = Names, y = Estimate)) +
  geom_pointrange(aes(ymin = Low, ymax = High)) +
  geom_hline(yintercept = 0, linetype = 2, color = "grey60") +
  ylab("Coefficient (Log Odds)") + xlab("") +
  theme_classic() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

# clear
rm(temp, ratio_data, logit_ratio, coef_plot_data, crossmove_data)


## 2: Analyses related to clique membership and appointment outcomes ======

favored_data <- read.csv("clique_favored_data.csv", header = T)

# == == == == === == === == == == == == === === == ==  == == == == == === == == ==
# == == == == Figure 5: Appointment outcomes across == == == == == == == == == === 
# == == == ==           clique and non-clique Officials == == == == === == == == =
# == == == == === == === == == == == == === === == ==  == == == == == === == == ==

figure5 <- favored_data

# Convert 1 to Clique and 0 to Non-clique to make visualization clear
figure5$Faction_Member <- ifelse(figure5$Faction_Member == 1, 
                                 "Clique", "Non_clique")
figure5$Faction_Member <- factor(figure5$Faction_Member)

# Create column to indicate � 1 and � 2 Bandwidth 
bandwidth <- append(rep("\u00B1 1 Bandwidth",
                        nrow(figure5[figure5$Year >= 60 & figure5$Year <= 61 , ])),
                    rep("\u00B1 2 Bandwidth",
                        nrow(figure5[figure5$Year >= 59 & figure5$Year <= 62 , ])))

# Collate figure5 and bandwidth for ggplot display
plot_data <- cbind(rbind(figure5[figure5$Year >= 60 & figure5$Year <= 61 , ],
                         figure5[figure5$Year >= 59 & figure5$Year <= 62 , ]), bandwidth)


ggplot(plot_data, 
       aes(x = promotion, fill = Faction_Member)) +
  geom_bar(stat="count",width=0.3,  position=position_dodge()) +
  facet_grid(.~bandwidth) +
  theme_minimal() +
  xlab("") + ylab("Count") +
  scale_fill_grey() +
  theme( legend.position = "bottom", legend.title = element_blank())

rm(figure5, plot_data, bandwidth)

# == == == == === == === == == == == == === === == ==  == == == == == === == == ==
# == == == ==                                                    = == === == == ==
# == == == == Table 3: Clique vs Non-clique Appointment Outcomes == == == == == == 
# == == == ==                                                    = == === == == ==
# == == == == === == === == == == == == === === == ==  == == == == == === == == ==

# Run models with 1 year bandwidth
favor_1 <- glm(favored ~ Faction_Member + Postcoup +
                 Faction_Member*Postcoup
               + Experience + Seq + Time_Past,
               data = favored_data[favored_data$Year >= 60 & favored_data$Year <= 61 , ], 
               family = "binomial")

disfavor_1 <- glm(disfavored ~ Faction_Member + Postcoup +
                    Faction_Member*Postcoup
                  + Experience + Seq + Time_Past,
                  data = favored_data[favored_data$Year >= 60 & favored_data$Year <= 61 , ], 
                  family = "binomial")

neutral_1 <- glm(neutral ~ Faction_Member + Postcoup +
                   Faction_Member*Postcoup
                 + Experience + Seq + Time_Past,
                 data = favored_data[favored_data$Year >= 60 & favored_data$Year <= 61 , ], 
                 family = "binomial")




# 1 Year bandwidth results
stargazer(favor_1, disfavor_1, neutral_1,
          covariate.labels=c("Clique Member", "Post Coup","Experience",
                             "Sequence", "Time Elapsed", "Clique Member*Post Coup"),
          style = "AJPS", type = "text")



rm(disfavor_1, favor_1, neutral_1)

# == == == == == === == == = = === == == == == == == == == == === == = == == == #
# === == == == == == == == == Optional: Run Models with == == === == === == === #
# == == == === == === == == =           2 Year bandwidth == == === === == == == #
# == == == == == === == == = = === == == == == == == == == == === == = == == == #
# Run same models with 2 year bandwidth
# favor_2 <- glm(favored ~ Faction_Member + Postcoup +
#                  Faction_Member*Postcoup
#                + Experience + Seq + Time_Past,
#                data = favored_data[favored_data$Year >= 59 & favored_data$Year <= 62 , ], 
#                family = "binomial")
# 
# disfavor_2 <- glm(disfavored ~ Faction_Member + Postcoup +
#                     Faction_Member*Postcoup
#                   + Experience + Seq + Time_Past,
#                   data = favored_data[favored_data$Year >= 59 & favored_data$Year <= 62 , ], 
#                   family = "binomial")
# 
# neutral_2 <- glm(neutral ~ Faction_Member + Postcoup +
#                    Faction_Member*Postcoup
#                  + Experience + Seq + Time_Past,
#                  data = favored_data[favored_data$Year >= 59 & favored_data$Year <= 62 , ], 
#                  family = "binomial")


# stargazer(favor_2, disfavor_2, neutral_2,
#           covariate.labels=c("Clique Member", "Post Coup","Experience",
#                              "Sequence", "Time Elapsed", "Clique Member*Post Coup"),
#           style = "AJPS", type = "text")




# == == == == == === == == = = === == == == == == == == == == === == = == == == #
# === == == == == == == == == Set up First difference simulation == == === == = #
# == == == === == === == == = to be used for Figure 6 and 7      === === == =  #
# == == == == == === == == = = === == == == == == == == == == === == = == == == #


first.diff.sim <- function(input_data, simulation_count, logit_outcome){
  
  # Objective:
  # 1. Clique membership is the 'treatment'
  # 2. Run a logit model to determine its correlation to the outcome of interest
  # 3. Perform the necessary simulations to obtain the FIRST DIFFERENCE
  #    between the predicted probability where clique membership is 1 vs 0
  
  
  #== == == === == == === == Step 1 == == === == == === == == === = =#
  set.seed(231)
  
  sim_mean <- NULL
  sim_sd <- NULL
  input_data <- input_data
  
  dv_set <- cbind(Intercept = 1, 
                  with(input_data, 
                       data.frame(Faction_Member = 0, Postcoup = 0,
                                  Experience = mean(Experience), 
                                  Seq = mean(Seq), Time_Past = mean(Time_Past),
                                  Faction_Member_Postcoup = 0)))
  
  dv_set <- rbind(dv_set, dv_set, dv_set, dv_set)
  dv_set[2, 2] <- 1 
  dv_set[3, 3] <- 1
  dv_set[4, c(2,7)] <- 1
  
  #== == == === == == === == Step 2 == == === == == === == == === = =#
  
  # Run the logit
  if(logit_outcome == "favored"){
    
    logit.m <- glm(favored ~ Faction_Member + Postcoup +
                     Faction_Member*Postcoup +
                     Experience + Seq + Time_Past,
                   data = input_data, 
                   family = "binomial")
    
  }else if(logit_outcome == "disfavored"){
    
    logit.m <- glm(disfavored ~ Faction_Member + Postcoup +
                     Faction_Member*Postcoup +
                     Experience + Seq + Time_Past,
                   data = input_data, 
                   family = "binomial")
    
  }else if(logit_outcome == "neutral"){
    
    logit.m <- glm(neutral ~ Faction_Member + Postcoup +
                     Faction_Member*Postcoup +
                     Experience + Seq + Time_Past,
                   data = input_data, 
                   family = "binomial")
  }
  
  
  # Simulate coefficients from the logit
  sim_logit <- mvrnorm(simulation_count, coef(logit.m), vcov(logit.m))
  
  
  # Matrix multiplication to perform linear combination 
  control_sim_prior <- sim_logit %*% t(dv_set[1, ]) # Clique = 0, Coup = 0
  clique_sim_prior <- sim_logit %*% t(dv_set[2, ])   # Clique = 1, Coup = 0
  
  control_sim_after <- sim_logit %*% t(dv_set[3, ]) # Clique = 0, Coup = 1
  clique_sim_after <- sim_logit %*% t(dv_set[4, ])   # Clique = 1, Coup = 1
  
  # Turn into probabilities
  control_sim_prior <- exp(control_sim_prior)/(1+exp(control_sim_prior)) 
  clique_sim_prior <- exp(clique_sim_prior)/(1+exp(clique_sim_prior)) 
  
  control_sim_after <- exp(control_sim_after)/(1+exp(control_sim_after)) 
  clique_sim_after <- exp(clique_sim_after)/(1+exp(clique_sim_after))
  
  
  # == == === == == Perform Difference in Difference == === == === == ==== = #
  
  control_control <- control_sim_after - control_sim_prior
  treat_treat <- clique_sim_after - clique_sim_prior
  diff_diff <- treat_treat - control_control
  
  treat_control_prior <- clique_sim_prior - control_sim_prior
  treat_control_after <- clique_sim_after - control_sim_after
  
  
  final_output <- as.data.frame(cbind(c("control_control","treat_treat", 
                                        "diff_diff","treat_control_prior", 
                                        "treat_control_after"), 
                                      rbind(cbind(mean(control_control), 
                                                  sd(control_control)),
                                            cbind(mean(treat_treat), 
                                                  sd(treat_treat)),
                                            cbind(mean(diff_diff), 
                                                  sd(diff_diff)),
                                            cbind(mean(treat_control_prior), 
                                                  sd(treat_control_prior)),
                                            cbind(mean(treat_control_after), 
                                                  sd(treat_control_after)))))
  
  
  
  
  colnames(final_output) <- c("group","mean", "sd")
  
  final_output$mean <- as.numeric(as.character(final_output$mean))
  final_output$sd <- as.numeric(as.character(final_output$sd))
  
  return(final_output)
  
}

#== === == === === === ==    First Difference         == == === == === === ==

# Obtain values for each outcome: Disfavored, Favored, Neutral

# 1 year bandwidth
dis <- first.diff.sim(favored_data[favored_data$Year >= 60 & favored_data$Year <= 61 , ],
                      1000, "disfavored")
fav <- first.diff.sim(favored_data[favored_data$Year >= 60 & favored_data$Year <= 61 , ],
                      1000, "favored")
neutral <- first.diff.sim(favored_data[favored_data$Year >= 60 & favored_data$Year <= 61 , ],
                          1000, "neutral")

# 2 year bandwidth
dis_62 <- first.diff.sim(favored_data[favored_data$Year >= 59 & favored_data$Year <= 62 , ],
                         1000, "disfavored")
fav_62 <- first.diff.sim(favored_data[favored_data$Year >= 59 & favored_data$Year <= 62 , ],
                         1000, "favored")
neutral_62 <- first.diff.sim(favored_data[favored_data$Year >= 59 & favored_data$Year <= 62 , ],
                             1000, "neutral")


#=== === === = == === == === === == === == === === === === === === 
#==                                                             ==
#= Figure 6: Pre-post coup difference among clique (left panel) ==
#==          and non-clique officials (right panel)             ==
#==          as generated from Table 3                          ==
#==                                                             ==
#== == === === === == === == === === === == === == == == === === =

outcome <-  c("Disfavored", "Disfavored",
              "Disfavored", "Disfavored",
              "Favored",  "Favored",
              "Favored",  "Favored",
              "Neutral", "Neutral", 
              "Neutral", "Neutral")

time <- c("\u00B1 1 Bandwidth",
          "\u00B1 1 Bandwidth",
          "\u00B1 2 Bandwidth",
          "\u00B1 2 Bandwidth", 
          "\u00B1 1 Bandwidth", 
          "\u00B1 1 Bandwidth",
          "\u00B1 2 Bandwidth", 
          "\u00B1 2 Bandwidth", 
          "\u00B1 1 Bandwidth",
          "\u00B1 1 Bandwidth",
          "\u00B1 2 Bandwidth", 
          "\u00B1 2 Bandwidth")



diff_one <- as.data.frame(cbind(time, outcome,
                                rbind(dis[1:2, ],  
                                      dis_62[1:2, ],
                                      fav[1:2, ],
                                      fav_62[1:2, ], 
                                      neutral[1:2, ], 
                                      neutral_62[1:2, ]
                                )))


diff_one[ , 3] <- factor(c("Non-clique (Post) - Non-clique (Pre)",
                           "Clique (Post) - Clique (Pre)"))



diff_one_1 <- diff_one[c(1:2, 5:6, 9:10) , ]
diff_one_2 <- diff_one[-c(1:2, 5:6, 9:10) , ]


# Plot 1 year bandwidth
ggplot(diff_one_1, aes(x = outcome, y = mean)) +
  geom_pointrange(aes(ymin = mean - 1.96*sd, ymax = mean + 1.96*sd), 
                  lwd = 1/2) +
  geom_linerange(aes(x = outcome, ymin = mean - 1.64*sd, ymax = mean + 1.64*sd),
                 lwd = 1) +
  geom_hline(yintercept = 0, linetype = 2, color = "grey60") +
  ylab("Difference in Predicted Probability") + xlab("") +
  facet_grid(.~group) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

# Optional: 2 year bandwidth
# ggplot(diff_one_2, aes(x = outcome, y = mean)) +
#   geom_pointrange(aes(ymin = mean - 1.96*sd, ymax = mean + 1.96*sd), 
#                   lwd = 1/2) +
#   geom_linerange(aes(x = outcome, ymin = mean - 1.64*sd, ymax = mean + 1.64*sd),
#                  lwd = 1) +
#   geom_hline(yintercept = 0, linetype = 2, color = "grey60") +
#   ylab("Difference in Predicted Probability") + xlab("") +
#   facet_grid(.~group) +
#   theme_bw() +
#   theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())


rm(outcome, time, diff_one, diff_one_1, diff_one_2)


#=== === === = == === == === === == === == === === === ===  === 
#==                                                          ==
#= Figure 7: First difference in the predicted probabilities ==
#==          for clique versus non-clique members            ==
#==          as generated from Table 3                       ==
#==                                                          ==
#=== === === = == === == === === == === == === === === === == = 


outcome <-  c("Disfavored", "Disfavored", "Disfavored", "Disfavored",
              "Favored", "Favored", "Favored", "Favored",
              "Neutral", "Neutral", "Neutral", "Neutral")

time <- c("\u00B1 1 Bandwidth", "\u00B1 1 Bandwidth", 
          "\u00B1 2 Bandwidth", "\u00B1 2 Bandwidth", 
          "\u00B1 1 Bandwidth", "\u00B1 1 Bandwidth", 
          "\u00B1 2 Bandwidth", "\u00B1 2 Bandwidth",
          "\u00B1 1 Bandwidth", "\u00B1 1 Bandwidth", 
          "\u00B1 2 Bandwidth", "\u00B1 2 Bandwidth")

model <- c("Pre Coup", "Post Coup",
           "Pre Coup", "Post Coup", 
           "Pre Coup", "Post Coup", 
           "Pre Coup", "Post Coup", 
           "Pre Coup", "Post Coup",
           "Pre Coup", "Post Coup")

temp <- as.data.frame(cbind(time, model, outcome,
                            rbind(dis[c(4,5), ],  
                                  dis_62[c(4,5), ],
                                  fav[c(4,5), ],
                                  fav_62[c(4,5), ], 
                                  neutral[c(4,5), ], 
                                  neutral_62[c(4,5), ]
                            )))


temp$model <- factor(temp$model, levels = c("Pre Coup", "Post Coup"))


# Reduce to just 1 Year time-bin
clique_non_diff_1 <- temp[temp$time == "\u00B1 1 Bandwidth", ]

# Expand to 2 Year time-bin
clique_non_diff_2 <- temp[temp$time == "\u00B1 2 Bandwidth", ]


# Plot 1 year bandwidth
ggplot(clique_non_diff_1, aes(x = outcome, y = mean)) +
  geom_pointrange(aes(ymin = mean - 1.96*sd, ymax = mean + 1.96*sd, color = model), 
                  lwd = 1/2, position = position_dodge(0.4)) +
  geom_linerange(aes(x = outcome, ymin = mean - 1.64*sd, ymax = mean + 1.64*sd, color = model),
                 lwd = 1, position = position_dodge(0.4)) +
  geom_hline(yintercept = 0, linetype = 2, color = "grey60") +
  ylab("Difference in Predicted Probability") + xlab("") +
  scale_color_manual(values=c('dark gray', 'black', 'dark gray', 'black'))+
  theme_classic() +
  theme(legend.title=element_blank(), legend.position = c(.518, .93),
        legend.direction = "horizontal", legend.text=element_text(size= 7))

# Optional: 2 year bandwidth
# ggplot(clique_non_diff_2, aes(x = outcome, y = mean)) +
#   geom_pointrange(aes(ymin = mean - 1.96*sd, ymax = mean + 1.96*sd, color = model), 
#                   lwd = 1/2, position = position_dodge(0.4)) +
#   geom_linerange(aes(x = outcome, ymin = mean - 1.64*sd, ymax = mean + 1.64*sd, color = model),
#                  lwd = 1, position = position_dodge(0.4)) +
#   geom_hline(yintercept = 0, linetype = 2, color = "grey60") +
#   ylab("Difference in Predicted Probability") + xlab("") +
#   scale_color_manual(values=c('dark gray', 'black', 'dark gray', 'black'))+
#   theme_classic() +
#   theme(legend.title=element_blank(), legend.position = c(.518, .93),
#         legend.direction = "horizontal", legend.text=element_text(size= 7))

# Optional: 1 and 2 year bandwidth side by side
# ggplot(temp, aes(x = outcome, y = mean)) +
#   geom_pointrange(aes(ymin = mean - 1.96*sd, ymax = mean + 1.96*sd, color = model), 
#                   lwd = 1/2, position = position_dodge(0.4)) +
#   geom_linerange(aes(x = outcome, ymin = mean - 1.64*sd, ymax = mean + 1.64*sd, color = model),
#                  lwd = 1, position = position_dodge(0.4)) +
#   geom_hline(yintercept = 0, linetype = 2, color = "grey60") +
#   ylab("Difference in Predicted Probability") + xlab("") +
#   scale_color_manual(values=c('dark gray', 'black', 'dark gray', 'black'))+
#   facet_grid(.~time) +
#   theme_bw() +
#   theme(legend.title=element_blank(), legend.position = c(.818, .93),
#         legend.direction = "horizontal", legend.text=element_text(size= 7))
